import numpy as np
from skimage import data,img_as_float,img_as_ubyte,exposure,io,color
from skimage.io import imread
from skimage.exposure import cumulative_distribution
from skimage.restoration import denoise_bilateral,denoise_nl_means,estimate_sigma
from skimage.measure import compare_psnr
from skimage.util import random_noise
from skimage.color import rgb2gray
from PIL import Image,ImageEnhance,ImageFilter
from scipy import ndimage,misc
import matplotlib.pylab as plt
import matplotlib

plt.rcParams['font.sans-serif']=['simhei']
plt.rcParams['font.family']='sans-serif'

# 4.1.1
# 实现输入图像的颜色通道直方图
def plot_image(image,title=''):
    plt.title(title,size=20),plt.imshow(image)
    plt.axis('off')
def plot_hist(r,g,b,title=''):
    r,g,b=img_as_ubyte(r),img_as_ubyte(g),img_as_ubyte(b)
    plt.hist(np.array(r).ravel(),bins=256,range=(0,256),color='r',alpha=0.5)
    plt.hist(np.array(g).ravel(),bins=256,range=(0,256),color='g',alpha=0.5)
    plt.hist(np.array(b).ravel(),bins=256,range=(0,256),color='b',alpha=0.5)
    plt.xlabel('像素值',size=20)
    plt.ylabel('频率',size=20)
    plt.title(title,size=20)
im=Image.open('parrot.png')
im_r,im_g,im_b=im.split()
plt.style.use('ggplot')
plt.figure(figsize=(15,5))
plt.subplot(121),plot_image(im,'原始图像')
plt.subplot(122),plot_hist(im_r,im_g,im_b,'RGB颜色通道直方图')
plt.show()

# 使用PIL图像模块的point（）函数进行对数变换，并转换作用于RGB图像，对不同色彩通道直方图的影响
im=im.point(lambda i:255*np.log(1+i/255))
im_r,im_g,im_b=im.split()
plt.style.use('ggplot')
plt.figure(figsize=(15,5))
plt.subplot(121),plot_image(im,'经对数变换后的图像')
plt.subplot(122),plot_hist(im_r,im_g,im_b,'RGB色彩通道对数变换直方图')
plt.show()


# 4.1.2
# 应用幂律变换对具有scikit—image的rgb进行变换，然后可视化变换对颜色通道直方图的影响
im=img_as_float(imread('eagle.png'))
gamma=5
im1=im**gamma
plt.style.use('ggplot')
plt.figure(figsize=(15,5))
plt.subplot(121),plot_hist(im[...,0],im[...,1],im[...,2],'RGB颜色通道（输入）直方图')
plt.subplot(122),plot_hist(im1[...,0],im[...,1],im[...,2],'RGB颜色通道（输出）直方图')
plt.show() 


# 4.1.3
# 将RGB图像划分为不同颜色通道，以可视化不同颜色通道像素的直方图
im=Image.open('cheetah.png')
im_r,im_g,im_b,_=im.split()
plt.style.use('ggplot')
plt.figure(figsize=(15,5))
plt.subplot(121)
plot_image(im)
plt.subplot(122)
plot_hist(im_r,im_g,im_b)
plt.show()

# PIL的point（）函数实现对比度拉伸
def contrast(c):
    return 0 if c < 70 else (255 if c > 150 else (255*c-22950)/48)

im1=im.point(contrast)
im_r,im_g,im_b,_=im1.split()
plt.style.use('ggplot')
plt.figure(figsize=(15,5))
plt.subplot(121)
plot_image(im1)
plt.subplot(122)
plot_hist(im_r,im_g,im_b)
plt.yscale('log',basey=10)
plt.show()

#使用PIL的ImageEnhance 模块
contrast=ImageEnhance.Contrast(im)
im1=np.reshape(np.array(contrast.enhance(2).getdata()).astype(np.uint8),(im.height,im.width,4))
plt.style.use('ggplot')
plt.figure(figsize=(15,5))
plt.subplot(121)
plot_image(im1)
plt.subplot(122)
plot_hist(im1[...,0],im1[...,1],im1[...,2])
plt.yscale('log',basey=10)
plt.show()


# 4.1.4
# 固定阈值的二值化
im=Image.open('cat1.png').convert('L')
plt.hist(np.array(im).ravel(),bins=256,range=(0,256),color='g')
plt.xlabel('Pixel values')
plt.ylabel('Frequency')
plt.title('Histogram of pixel values')
plt.show()
plt.figure(figsize=(12,18))
plt.gray()
plt.subplot(221),plot_image(im,'original image')
plt.axis('off')
th=[0,50,100,150,200]
for i in range(2,5):
    im1=im.point(lambda x:x>th[i],'1')
    plt.subplot(2,2,i),plot_image(im1,'binary image with threshold='+str(th[i]))
plt.show()

# 半色调二值化
im=Image.open('C:/Users/TCJ最帅/Desktop/tcj/images/cat1.png').convert('L')
im=Image.fromarray(np.clip(im+np.random.randint(-128,128,(im.height,im.width)),0,255).astype(np.uint8))
plt.figure(figsize=(12,18))
plt.subplot(221),plot_image(im,'带有噪声的原始图像')
th=[0,50,100,150,200]
for i in range(2,5):
    im1=im.point(lambda x:x>th[i],'1')
    plt.subplot(2,2,i),plot_image(im1,'阈值为'+str(th[i])+'的二值图像')
plt.show()


# 4.2.1
# 使用曝光模块的equalize_hist（）函数对scikit-image进行直方图均衡化
img=rgb2gray(imread('goal1.png'))
img_eq=exposure.equalize_hist(img)
img_adapteq=exposure.equalize_adapthist(img,clip_limit=0.03)
plt.gray()
images=[img,img_eq,img_adapteq]
titles=['原始输入图像','直方图均衡化后的图像','自适应直方图均衡化后的图像']
for i in range(3):
    plt.figure(figsize=(20,10)),plot_image(images[i],titles[i])
plt.figure(figsize=(15,5))
for i in range(3):
    plt.subplot(1,3,i+1)
    plt.hist(images[i].ravel(),color='g')
    plt.title(titles[i],size=15)
plt.show()    

# 两种不同的直方图处理技术，基于scikit-image的对比度拉伸和直方图均衡化得到的图像增强进行比较
matplotlib.rcParams['font.size']=8
def plot_image_and_hist(image,axes,bins=256):
    image=img_as_float(image)
    axes_image,axes_hist=axes
    axes_cdf=axes_hist.twinx()
    axes_image.imshow(image,cmap=plt.cm.gray)
    axes_image.set_axis_off()
    axes_hist.hist(image.ravel(),bins=bins,histtype='step',color='black')
    axes_hist.set_xlim(0,1)
    axes_hist.set_xlabel('Pixel intensity',size=15)
    axes_hist.ticklabel_format(axis='y',style='scientific',scilimits=(0,0))
    axes_hist.set_yticks([])
    image_cdf,bins=exposure.cumulative_distribution(image,bins)
    axes_cdf.plot(bins,image_cdf,'r')
    axes_cdf.set_yticks([])
    return axes_image,axes_hist,axes_cdf
    
im_feng=io.imread('cheetah.png')
im_f_rescale=exposure.rescale_intensity(im_feng,in_range=(0,100),out_range=(0,255))
im_f_eq=exposure.equalize_hist(im_feng)
im_adapteq=exposure.equalize_adapthist(im_feng,clip_limit=0.03) 

fig=plt.figure(figsize=(15,7))
axes=np.zeros((2,4),dtype=np.object)
axes[0,0]=fig.add_subplot(2,4,1)
for i in range(1,4):
    axes[0,i]=fig.add_subplot(2,4,1+i,sharex=axes[0,0],sharey=axes[0,0])
for i in range(0,4):
    axes[1,i]=fig.add_subplot(2,4,5+i)
axes_image,axes_hist,axes_cdf=plot_image_and_hist(im_feng,axes[:,0])
axes_image.set_title('低对比度图像',size=20)
y_min,y_max=axes_hist.get_ylim()
axes_hist.set_ylabel('像素总数量',size=20)
axes_hist.set_yticks(np.linspace(0,y_max,5))
axes_image,axes_hist,axes_cdf=plot_image_and_hist(im_f_rescale,axes[:,1])
axes_image.set_title('对比度拉伸',size=20)
axes_image,axes_hist,axes_cdf=plot_image_and_hist(im_f_eq,axes[:,2])
axes_image.set_title('直方图均衡化',size=20)
axes_image,axes_hist,axes_cdf=plot_image_and_hist(im_adapteq,axes[:,3])
axes_image.set_title('自适应均衡化',size=20)
axes_cdf.set_yticks(np.linspace(0,1,5))
fig.tight_layout()
plt.show()


# 4.2.2
# 直方图匹配
def cdf(im):
    c,b=cumulative_distribution(im)
    c=np.insert(c,0,[0]*b[0])
    c=np.append(c,[1]*(255-b[-1]))
    return c
def hist_matching(c,c_t,im):
    pixels =np.arange(256)
    new_pixels= np.interp(c,c_t,pixels)
    im =(np.reshape(new_pixels[im.ravel()],im.shape)).astype(np.uint8)
    return im

plt.gray()
im =(rgb2gray (imread ('beans_g.png'))*255).astype(np.uint8)
im_t =(rgb2gray (imread('lena_g.png'))*255).astype(np.uint8)

plt.figure(figsize=(20,12))
plt.subplot (2,3,1),plot_image(im,'输入图像')
plt.subplot(2,3,2),plot_image (im_t,'模板图像')
c=cdf(im)
c_t =cdf(im_t)
plt.subplot(2,3,3)
p= np.arange(256)
plt.plot (p,c,'r.-',label='输入图像')
plt.plot (p,c_t,'b.-',label='模板图像')
plt.legend(prop={'size':15})
plt.title('累积积分函数',size=20)
im= hist_matching(c,c_t,im)
plt.subplot(2,3,4),plot_image(im,'直方图匹配的输出图像')
c1 =cdf(im)
plt.subplot(2,3,5)
plt.plot(np.arange(256),c,'r.-',label='输入图像')
plt.plot(np.arange(256),c_t,'b.-',label='模板图像')
plt.plot(np.arange(256),c1,'g.-',label='输出图像')
plt.legend(prop={'size':15})
plt.title('累积积分函数',size=20)
plt.show()


# 4.3.1
# 基于ImageFilter 平滑
i=1
plt.figure(figsize=(10,25))
for prop_noise in np.linspace(0.05,0.3,3):
    im= Image.open('mandrill.jpg')
    n =int(im.width*im.height*prop_noise)
    x,y =np.random.randint (0,im.width,n),np.random.randint (0,im.height,n)
    for (x,y) in zip(x,y):
        im.putpixel((x,y),((0,0,0) if np.random.rand()<0.5 else (255,255,255)))
    im.save('mandrill_spnoise_'+str(prop_noise)+'.jpg')
    plt.subplot (6,2,i),plot_image(im,'掺杂了'+str(int(100*prop_noise))+'% 噪声的原始图像')
    i+=1
    im1= im.filter(ImageFilter.BLUR)
    plt.subplot (6,2,i),plot_image (im1,'模糊图像')
    i+=1
plt.show()

# 基于盒模糊核均质化平滑
im =Image.open ('mandrill_spnoise_0.3.jpg')
plt.figure(figsize=(20,7))
plt.subplot (1,3,1)
plt.imshow(im)
plt.title ('原始图像',size=30)
plt.axis ('off')
for n in [3,5]:
    box_blur_kernel= np.reshape(np.ones (n*n),(n,n))/(n*n)
    im1= im.filter(ImageFilter.Kernel((n,n),box_blur_kernel.flatten()))
    plt.subplot (1,3,(2 if n==3 else 3))
    plot_image(im1,'核大小为'+str(n)+'x'+str(n)+'的模糊化')
plt.suptitle('利用不同核大小的PIL均值滤波（盒模糊）',size=30)
plt.show()

# 基于高斯模糊滤波器平滑
im =Image.open('mandrill_spnoise_0.05.jpg')
plt.figure(figsize=(20,6))
i=1

for radius in range(1,4):
    im1= im.filter(ImageFilter.GaussianBlur(radius))
    plt.subplot(1,3,i),plot_image(im1,'半径为'+str (round (radius,2)))
    i+=1
plt.suptitle('不同半径的PIL高斯模糊',size=20)
plt.show()


# 4.3.2
#使用SciPy的ndimage模块函数对图像进行线性滤波
im =imread('mandrill_spnoise_0.3.jpg')
k =7
im_box= ndimage.uniform_filter(im,size=(k,k,1))
s =2 
t =(((k -1)/2)-0.5)/s 
im_gaussian= ndimage.gaussian_filter(im,sigma=(s,s,0),truncate=t)
fig=plt.figure(figsize=(30,10))
plt.subplot(131),plot_image(im,'原始图像')
plt.subplot(132),plot_image(im_box,'带有盒式滤波器（的模糊化）')
plt.subplot(133),plot_image(im_gaussian,'带有高斯滤波器（的模糊化）')
plt.show()


# 4.4.1
# 中值滤波器
i=1
plt.figure(figsize=(25,35))
for prop_noise in np.linspace(0.05,0.3,3):
    im =Image.open('mandrill.jpg')
    n =int (im.width *im.height *prop_noise)
    x,y =np.random.randint (0,im.width,n),np.random.randint (0,im.height,n)
    for (x,y) in zip(x,y):
        im.putpixel((x,y),((0,0,0) if np.random.rand()<0.5 else (255,255,255)))
    im.save('mandrill_spnoise_'+str(prop_noise)+'.jpg')
    plt.subplot(6,4,i)
    plot_image(im,'掺杂了 '+str(int (100*prop_noise))+'%噪声的原始图像')
    i+=1
    for sz in [3,7,11]:
        im1= im.filter(ImageFilter.MedianFilter(size=sz))
        plt.subplot(6,4,i),plot_image(im1,'输出(中值滤波器大小为'+str(sz)+')')
        i+=1
plt.show()

# 最大值滤波器核最小值滤波器
im =Image.open('mandrill_spnoise_0.3.jpg')
plt.subplot(1,3,1)
plot_image(im,'掺杂了30%噪声的原始图像')
im1= im.filter(ImageFilter.MaxFilter(size=sz))
plt.subplot(1,3,2),plot_image(im1,'Output (Max Filter size='+str(sz)+')')
im1=im1.filter(ImageFilter.MinFilter(size=sz))
plt.subplot(1,3,3),plot_image(im1,'Output (Min Filter size='+str(sz)+')',size=15)
plt.show()


# 4.4.2
# 双边滤波器
im=color.rgb2gray(img_as_float(io.imread('hill.png')))
sigma=.0155
noisy=random_noise(im,var=sigma**2)
plt.imshow(noisy)
plt.show()

# 使用双边滤波器对噪声图像去噪
plt.figure(figsize=(20,15))
i=1
for sigma_sp in [5,10,20]:
    for sigma_col in [0.1,0.25,5]:
        plt.subplot (3,3,i)
        plt.imshow(denoise_bilateral (noisy,sigma_color=sigma_col,sigma_spatial=sigma_sp,multichannel=False))
        plt.title(r'$\sigma_r=$'+str(sigma_col)+r',$\sigma_s=$'+str(sigma_sp),size=20)
        i+=1
plt.show()

# 非局部均值滤波器
def plot_image_axes (image,axes,title):
    axes.imshow (image)
    axes.axis ('off')
    axes.set_title(title,size=20)
parrot=img_as_float (imread('parrot.png'))
sigma =0.25
noisy=parrot +sigma *np.random.standard_normal(parrot.shape)
noisy=np.clip(noisy,0,1)
sigma_est =np.mean(estimate_sigma (noisy,multichannel=True))
print("estimated noise standard deviation ={}".format(sigma_est))
patch_kw=dict(patch_size=5,patch_distance=6,multichannel=True)
denoise=denoise_nl_means(noisy,h=1.15*sigma_est,fast_mode=False,**patch_kw)
denoise_fast=denoise_nl_means(noisy,h=0.8*sigma_est,fast_mode=True,**patch_kw)
fig,axes=plt.subplots(nrows=2,ncols=2,figsize=(15,12),sharex=True,sharey=True)
plot_image_axes(noisy,axes[0,0],'噪声图像')
plot_image_axes(denoise,axes[0,1],'非局部均值（慢）')
plot_image_axes(parrot,axes[1,0],'原始图像（无噪声）')
plot_image_axes(denoise_fast,axes[1,1],'非局部均值（快）')
fig.tight_layout()
psnr_noisy=compare_psnr(parrot,noisy)
psnr=compare_psnr(parrot,denoise.astype(np.float64))
psnr_fast=compare_psnr(parrot,denoise_fast.astype(np.float64))
print("PSNR (noisy)={:0.2f}".format(psnr_noisy))
print ("PSNR (slow)={:0.2f}".format(psnr))
print ("PSNR (fast)={:0.2f}".format(psnr_fast))
plt.show()


# 4.4.3
# Scipy ndimage 平滑
lena=imread('lena.jpg')
noise=np.random.random(lena.shape)
lena[noise >0.9]=255
lena[noise <0.1]=0
plot_image (lena,'noisy image')
plt.show ()
fig=plt.figure(figsize=(20,15))
i=1
for p in range(25,100,25):
    for k in range(5,25,5):
        plt.subplot (3,4,i)
        filtered=ndimage.percentile_filter (lena,percentile=p,size=(k,k,1))
        plot_image(filtered,str(p)+'%核尺寸为'+str(k)+'x'+str(k))
        i+=1
plt.show ()